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Abstract 

In this article, we prove that Dirac brackets for Hamiltonian and non-Hamiltonian 
constrained systems can be derived recursively. We then study the applicability of 
that formulation in analysis of some interesting physical models. Particular attention 
is paid to feasibility of implementation code for Dirac brackets in Computer Algebra 
System and analytical techniques for inversion of triangular matrices. 
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1 Introduction 



The fundamental notion in the Hamiltonian formulation of classical dynamics 
of particles and fields is the canonical Poisson bracket defined over the space 
of all differentiable functions of the phase space (of even dimension), such 
that: for each two phase space functions f{q,p) and g{q,p) where {q,p) = 
{qi, . . . , g„,pi, . . . ,pn) denote generalized positions and momenta respectively. 
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This bracket is linear in each argument, skew-symmetric: {/,(?} = — {(?,/}, 
satisfies Leibniz identity: {f,g ■ h} = {f,g} ■ h + g ■ {f,h}, Jacobi iden- 
tity: {f,{g,h}} + {g,{hj}} + {h,{f,g}} = and is non-degenerate, i.e. 
if {f,g} = for all g, then / = const. This canonical Poisson bracket equips 
the phase space with a symplectic structure [I] . The Hamiltonian dynamics is 
then determined by defining the proper Hamiltonian function Ti. The evolution 
equation for any phase space function f{q,p) reads then: %: = % + {f^'H}. 

In applications one often encounters a situation when the phase space dy- 
namics is subject to certain external restricting conditions on the phase space 
variables called constraints. Often the constraints can be written in terms of 
some phase space functions (pi{q,p) = 0, and we will restrict our analysis to 
these cases only. The Hamiltonian formalism for such constrained systems re- 
quires modifications. These modifications have been first suggested by Dirac 
[2], and a brief account of the Dirac theory follows. 

Let (pi (with i = 1, . . . , L) denote all constraints for our Hamiltonian system. 
Those constraints can be divided into two classes by analyzing the Lx L skew- 
symmetric matrix of their mutual Poisson brackets Aij = (pj}. Since A is 
skew-symmetric, its rank K must be even. We assume that after relabeling of 
the (pi and/or redefining the constraints by taking their linear combinations 
(known as the Dirac separating constraints algorithm), the top left K x K 
submatrix of A, which we denote by W, is regular. The constraint functions 
(pK+i, ■ ■ ■ ,(pL are then called first class constraints, and are associated with lo- 
cal gauge symmetries [2], while (pi, . . . ,(pK are called second-class. In this work 
we will consider second-class constraints only, and for them we can introduce 
the Dirac bracket (DB)[2], of two phase space functions f,g: 



In the modern language of symplectic geometry, constrained Hamiltonian dy- 
namics can be represented by a triplet (M, A^, uj) where (M, u) is a symplectic 
manifold, namely Phase space, and is a constraint submanifold of M. The 
DB f ll.2p is the Poisson bracket on a symplectic submanifold A^' C A^, called 
second-class constraint manifold [HEIIIS]. 

Symplectic structure requires even dimensional manifolds and non-degenerate 
Poisson structure. Both these assumptions seem too restrictive and not al- 
ways applicable. With the appearance of non-canonical Poisson structure (PS) 
in rigid body dynamics, theory of magnetism, infinite dimensional PS in 
magneto-hydrodynamics, etc. and issues of geometric quantization, systematic 
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studies of the general Poisson bracket (PB) which is a Lie bracket satisfying 
the Leibniz identity, has become important. 

The fundamental geometric object in the description of any generalized Hamil- 
tonian dynamics is a Poisson manifold. Geometrically, Poisson manifold is a 
manifold endowed with a bivector field tt satisfying [vr, vr] = 0, where [■, ■] 
denotes the Schouten bracket [6] on multivector fields. Algebraically, M is a 
Poisson manifold if there is a Poisson bracket on the space of smooth functions 
defined on M. The Poisson bracket {■, ■} and the bivector field vr determine 
each other [5.T] by the formula {/,(?} = TT{df,dg). Both the geometric and 
algebraic characterization of Poisson manifolds are used in the literature. 

In the analysis of the constrained systems dynamics it is of predominant im- 
portance to formulate it as a usual Poisson structure on a submanifold of a 
non- constrained system's Poisson manifold. The conditions under which the 
Poisson structure on a submanifold is achievable was investigated in [8]l9] and 
the geometric derivation of the DB formula flL2l) via a procedure called geo- 
metric reduction of Poisson tensor was known [lOj. 

In many of the important physical applications the systems described are not 
purely Hamiltonian but also dissipative. The description of such combined 
dissipative-hamiltonian dynamics can be formulated in various ways, however 
one of them seems to be particularly elegant and allows to incorporate in it 
many methods developed in purely symplectic dynamics. This method was 
introduced first in the phase transformation kinetics in [TT] and then indepen- 
dently in p^][T3| and called metriplectic. The main point in metriplectic for- 
mulation [13] is that a mixed bracket obtained by adding a symmetric bracket 
to the Poisson bracket can successfully be used for description of dissipative 
systems. 

In the metriplectic framework, the underlying structure of a dissipative system 
consists of a Poisson and a symmetric bracket [13] , and the obvious generaliza- 
tion of this construction for constrained dissipative system (CDS) must consist 
of two DB ^14j: the usual skew-symmetric DB and the symmetric DB, which 
describe the Hamiltonian and dissipative part respectively. In [14] we have 
assumed that CDS be geometrically represented by a triplet {M,N,u — g), 
here is a submanifold of the symplectic manifold (M, uj) and g is a. covariant 
semimetric tensor. Generalized result can be easily obtained by replacing the 
symplectic 2-form u; by a contravariant Poisson tensor vr, and the covariant 
metric (0, 2) tensor g hj a. contravariant (semi/pseudo)-metric (2, 0) tensor G. 

The aim of the article is to give a formal (algebraic) proof of the recursiveness 
of symmetric and skew- symmetric DB. For the latter, this property probably 
has been known for years in practical calculation, but none algebraic proof 
seems to be available in the literature. The proof given in this paper is, to the 
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best of our knowledge, the first one. 

The paper is organized as follows. Section 2 resumes a construction leading 
to the DB-like formula in the general case and conditions of submanifold pos- 
sessing Poisson structure in the form of the DB. Section 3 presents rigorous 
proof for the recursiveness of symmetric and skew-symmetric DB. Section 
4 illustrates the constrained metriplectic formalism on two examples, using 
the computer algebra package Mathematica. Appendix A shows that sym- 
bolic/analytical difficulties appeared in the Dirac approach are unavoidable 
and that they also appear in the Lagrangian approach. Appendix B contains 
Dirac and LMM description for A^-pendulum, which serves as our numerical 
case study. Appendix C contains some techniques for analytical inversion of 
symmetric tridiagonal matrices, which we worked out in 2004. 

In this article, we denote a symmetric, skew-symmetric and general bracket 
by < ■, ■ >, {■, ■} and ri{-, ■) respectively. 



2 Geometric interpretation on Dirac-like brackets 

We begin by showing how an arbitrary i^'-tensor defined on a manifold M can 
be reduced in the (almost) Dirac sense to any submanifold of M, regardless 
of this tensor degeneracy. 

Conventionally we will denote the dual spaces to E , F, etc. and similar dual 
map to /, etc. by superscript asterisk, e.g E*, F* and /*, the annihilatoiH by 
superscript zero, e.g F^ and V^. Furthermore, denoting annihilation between 
elements of E and E* by (■ | ■), each bivector vr G A^E defines the map vr" : 
E* ^ E hj (7r^(C) \ ri) = TX^Q^rj) for t] E E*. The term almost Poisson 
structure means that this structure is bilinear and skew-symmetric, but does 
not necessarily satisfy the Jacobi identity. 

Let E he a linear space, F be its linear subspace and let E he a direct sum 
E = F(BV. This direct sum determines uniquely the projection p : E F and 
induces a splitting in its dual space E* = F*(BV* with F* = V'^ and V* = 
The direct decomposition on E* determines uniquely the map p* : F* —>■ E* 
which is the dual map of p. 

Definition 1 Each multilinear map K : (E*)^ R induces multilinear map 

Kp : (F*)'' ^ Rhy 

Kpiai, ■■■,ak)= K{p*{ai), ■ ■ ■,p*{ak)) . (2.1) 
^ The annihilator oi F C E is F'^ = {(j) e E* such that (/)(F) = 0} 
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We will call Kp an almost Dirac reduction of on F with respect to the 
direct sum E = F ®V (or with respect to the projection p). 

If K is symmetric or skew-symmetric then Kp has the same properties, but 
Kp may not inherit other algebraic properties of K. In particular, if K is 
non-negative, i.e. K{a, a, ■ ■ ■ , a) > 0, Va G E* , then Kp also is non-negative, 
but if i^' = TT is Poissonian, the bracket defined hy ttp may not satisfy the 
Jacobi identity. Thus, in general ttp is only almost Poissonian. The sufficient 
condition for T^p to be Poissonian is: 

Proposition 1 Suppose that n is a Poisson tensor in E such that F fl 7r^(F°) = 
{0}. Then np is a Poisson tensor. 

Proof. Since F fl 7r'*(F°) = {0}, one has 7r'*(C) G V for every ( G F°. Hence 
^(C) v) = (^''(C) I ^) = for every ( G F", rj G V'^. This orthogonality condition 
implies that n G A^E decomposes as tt = np + ny where Tip G A^F and 
TTy G A'^V. The identity [vr, vr] = implies that [np, np] = [up + tt , tc p — it] = 
— [27Tp + TTv, VTy]. Because of [up, np] G A^F C A'^EAF and —[2TTp + nv, vry] G 
A^-E A V, they must both be zeros. Therefore [np, up] = which means that 
TTp is Poissonian. 

Proposition ([T]) leads to the following concept of Dirac subspace: A linear 
subspace F of a Poisson space {E, n) is called a Dirac subspace if Fn7r''(F'') = 
{0}. Furthermore, since any projection p : E F defines a split E = F (BVp 
where Vp = [Id — p)F and the condition F fl tt^{F^) = {0} is equivalent to 
71^ (F^) C Vp, this condition suggests to introduce a concept of Dirac projection: 
A linear map p : F — > F is called a Dirac projection if p(7r'*(F'')) = 0. 

We now consider the non-linear case. Let us consider a smooth finite dimen- 
sional manifold M, a submanifold N (Z M and a regular distribution V on M 
(that is a smooth family of the subspaces of the tangent spaces, Vx C T^M) 
such that TxM = TxN (BVx for every x in A^. Thus V is complementary of TN 
in TM. 

For any k one-forms ai, ■ ■ ■ , ak, the reduction of {k, 0)-tensor field AT on A^ is 
defined by: 

Km (ai, - ■ ■ ,ak) = Kip*{ai), - ■ ■ ,p*iak)) . (2.2) 

We call the tensor field the almost Dirac reduction of K with respect to 
the submanifold A^ and the direct decomposition T^M = TN © V. 

Applying proposition ([1]) one gets the following 

Proposition 2 Let N be a submanifold of a Poisson manifold (M, vr). Sup- 
pose that TN n 7r''(TA^°) = {0} and tt^ is smooth. Then, ttn is a Poisson 



5 



tensor on N. 



Thus we obtain a sufficient condition for constructing Poisson structure on a 
submanifold. Furthermore, proposition ([2]) leads to the following concept of 
Dirac submanifold. 

Definition 2 A submanifold N of the Poisson manifold (M, it) is called a 
Dirac submanifold ifTNn tt^{TN^) = {0} and induced tensor is smooth. 

Note that the concept of Dirac submanifold (def. [2]) is less restrictive than 
the one introduced by Xu [8J and the other mentioned therein: A submanifold 
N is -called by Xu- a Dirac submanifold of the Poisson manifold (M, 11) if 
there exists a bundle V such that T^M = TN © V and V is a coisotropic 
submanifold of TM. 

For applications, the most important case of this geometric procedure is when 
K = IT ± G, where 7i,G are Poisson and pseudo/semi- metric tensor, respec- 
tively. 



3 Algebraic formulas for computing Dirac brackets 



3.1 Pfaffians and the Tanner's identities 



For any function of two arguments F defined on the set of generators of the 
commutative algebra A, we introduce the notation 



F[xi ■ ■ ■Xn,yi - ■ -iin] = det(F[a;i, yj]) 



F[xi,yi] ■ ■ ■ F[xi,yn] 
F[xn,yi] ■ ■ ■ F[xn,yn] 



(3.1) 



We will use the following identities: 



F[a, p] F[axz, fyt] = F[ax, py] F[az, pt] - F[ax, pt] F[az, Py] (3.2) 
F[a, P] F[axuv, Pyst] = F[ax, Py] F[auv, Pst] — F[ax, Ps] F[auv, Pyt] 

+ F[ax,pt]F[auv,pys], (3.3) 

which are a special case of the Tanner identity [T5p6] : and they also are 
known as theorems on bordered determinants [17], pages 46-50. Assuming 
F[u,v] = T]{u,v) for u,v from a commutative algebra with the bracket rj, we 
have 
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^(01,6) 



(3.4) 



3.2 Determinant and recursive formulas 



Let (JF, •) be a commutative algebra with the bracket ri : T ^ T T and 
{(/>i}"=i be a set of elements from T . Suppose the square matrix W = {Wij) 
with Wij = ri{(f)i, (pj) is invertible, and let us denote its inverse matrix by 
C = [Cij]. The original DB formula follows: 



The new bracket (13.51) is bilinear and it inherits algebraic properties from the 
original bracket rj. It is easy to check that V/ G J-", 'r]D{4'i, f) = 0, which means 
that all elements (pi are in the algebra center (called Casimir's elements) of the 
algebra (JF, rjr))- For skew-symmetric algebras the number of fixed elements (pj 
must be even, because the skew-symmetric matrix W with odd rank always 
is singular. Indeed, denoting det W by \W\, for skew-symmetric matrix W we 
have \W\ = \W^\ = (-l)"|Vr|. 

Let A = (aij) be a matrix, then the matrix obtained from A after deleting i— th 
row and j— th column will be denoted by A^^'^\ Recall the Laplace expansion 
formula which states that det A = \A\ = !]-,•(— l)*'''-'^^^!^^*'-'''! for any square 
matrix A. Now we can easily prove the following determinant formula for the 



Proposition 3 fl4^ Supposing the matrix W{(pi, . . . ,(pn) is invertible, the 
following identity holds 



n 



VD{f,9) = v{f,9) - v{f,<Pi)Cijri{(pj,g) , \/f,geJ^. 



(3.5) 



DB. 
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Wf,g eJ^ : VD{f,g) 



r/(0i,0i) ■ ■ ■ r/(0i,0„) v{(Pi,9) 



^(01, 0l) ■ • ■ Vih,<Pr. 



Rewriting (13.61) in the notation (13. ip we get 



Wf,g eJ^ : r]D{f,g) 



\W 



\W\ 



(3.6) 



(3.7) 



where \ W\ = F[ 



and \WfJ = F[(f)i ■ ■ ■ ■ ■ ■ (p^g]. 



Proof. Apply twice the Laplace formula to the last column and row of the 

matrix Wf^g. 

A. Symmetric case: 

Now let (JF, ■) be a commutative algebra with the bracket < ■, ■ > and {</>j}j=i, 
be a set of elements from JF. We define inductively a family of brackets 



<f,g>'^'^ = <f,g>, 



< >^^^< A+i,g > 

<0fc+i,0fc+i>W 



(fe) 



(3. 



Denote the Dirac bracket determined by k constraints (pa with a = 1, - ■ ■ ,k, 
< f,9 >d\ thus 



< f,9 >V=< f,g>-T. <f,<Pa> Clt> < <P,,g > 



(fc) 



(3.9) 



a, 6=1 



where C^'^^ is the inverse matrix oikxk matrix W^''^ 



We prove the following theorem 



8 



Theorem 1 (Recursive general brackets) Assume that the family of brack- 
ets ( 13. 8p is well-defined. Then V/, g & T and 1 < m < n: 

<f,g>^^^ = <f,g>t^ . (3.10) 



Proof. We prove the formula fl3.10p by induction with m. For m = 1, (13.101) 
is obviously true. Suppose that it is true for m = k, thus 

yf,g: (7 (3.11) 

we shall prove that it remains true for m = k + 1. The proof is based on the 
Tanner identity (13.21) and the proposition [31 

First, let a = 0i02 ■ ■ -(pk, using formula (13. 7p in the proposition |3] we have 



^ (fc+i)_ <Pk+if, a (pk+ig] ,oioN 
h [a (pk+i,a(pk+i\ 

Multiplying r.h.s. of (IXT^ by 1 = f}^ and using (13.21) we get 



^ ^ ^ ^ (fc+i) ^ F[af, ag] F[af, a 0fc+i] F[a 0fc+i, ag] 

Using formula fl3.7p again, we show that: the first term in the r.h.s. of eq. 03.131) 
is equal < f,g and also equal < f,g >^''^ by induction assumption (13. lip . 
Applying similar argument for the second term in the r.h.s. of eq. fl3.13p . we 
obtain 



] =< /, 0fc+i >^ ' , ] =< g >^ ' and 

— =< (pk+i,(pk+i >^ ' ■ 

r [a, a\ 

In summary, the r.h.s. of eq. 03.131) is equal 



< /, g >(^) _lA^±tl^!!l5i^if^ . (3.14) 

It implies that r.h.s. of eq. 03.130 is equal < f,g >(^+^) which ends the proof. ^ 
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To apply theorem [T] we need an existence of the family of brackets (13.81) . This 
condition requires the invertibility of < 4>i+i, (pi+i >^*^ for all i with 1 < i < n, 
and therefore it is equivalent to the regularity (or non- degeneracy) of all main 
minors of W. This condition may seem to be too restrictive, however by making 
new constraints from linear combinations of old constraints, we can go beyond 
this restriction The following simple example illustrates the procedure. 

Example 3.1 Let x = (xi, X2, ■ ■ ■ , x„) g -R", 



< Xi,Xi >=< X2, X2 >= 0, < Xi, X2 >=< X2, Xi >= a{x) , 
other brackets are whatever, and the constraints are 0i = xi = 0, 4)2 = X2 = 0. 



In the standard approach, after calculating the constraint matrix W = a{x) 
and its inverse, we easily get the Dirac bracket 

< f,9 >D=< f,g > — ^(< f,xi >< X2,g > + < f,X2 >< xi,g >) . 
In this case, direct recursive scheme is inapplicable because of 

<(/>!, 01 > = = <4>2,4>2>, 

but by introducing new (equivalent) constraints ui = xi + X2 = and U2 
Xi — X2 = 0, the recursive scheme may apply as below. 

In the first step, we have 
<f,g >« = < > 



1 

1 



(1) = < f > < >< ui,g > 

< Ml, Ml > 



Since < Ui,U2 >= we get < f,U2 >*^^^=< f,U2 >, < U2,g >^^^=< U2,g > 
and < U2,U2 >*^"'^-*=< ^2,^2 >• Hence, 



< U2,U2 >(^) 
<f,Ui><Ui,g> <f,U2><U2,g> 



= < f,g > 

< Ml, Ml > < M2, M2 > 

Finally, express it in terms of the original constraints 



10 



<f,9 >^^' = < f,9> — r^(< f,xi><X2,g > + < f,X2>< xi, g >) . 



We can use theorem [T] to prove that symmetric DB inherits non-negativity 
from a semimetric bracket. Precisely, 

Proposition 4 Suppose JF be an algebra of real functions with semimetric 
bracket <-,■>, i.e. < /, / > is a non-negative function for every function 
/ G JF. Let {(j)k}k=i be a set of elements from such that . . . , 0„) 

is invertible. Then the Dirac bracket < -,■ >d with respect to {(f)k}k=i, is 
semimetric. 

Proof. Since the recursion property of symmetric DB in theorem [T], it is 
enough to prove < /, / >*^^'' is a non-negative function. Indeed, for every real 
number A, one has 

< < / - X(jyiJ- Hi >=< fj> -2A </,</.!> +X' < <PiAi > , 

which implies that the discriminant A = [< /, 0i >]^— < f,f >< <t>i,<t>i >< 
0. Thus, 



< 01,01 > 

B. Skew-symmetric case: 

Now let {J^, ■) be a commutative algebra with a skew-symmetric bracket {■, ■} 
and {0fc}i=i5 be a set of elements from JF. We define inductively a family of 
brackets 



{/,5}^°^={/,5}, 

{/ ^jC^+i) = {/ - i/' 02fc+2}^^'H'/'2fc+i,g}^^^ - {/, 4>2k+iY^K4>2k+2.gY^^ 

(3.15) 

We prove that (13.151) are identical with the Dirac brackets. 

Theorem 2 (Recursive skew-symmetric brackets) Suppose that the fam- 
ily of bracket recursively defined by (13.151) is well-defined. Then V/, g ^ T and 

\ < m <n: 



(3.16) 
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where r.h.s. is the Dirac bracket with respect to 2m constraints 



2m 



{/, gVr = if, 9}-E {/> <l>a} 9} . 

a,b=l 



In the above C*^^*"-* in the inverse of the 2m x 2m matrix VT^^™-) 



{01, 0l} ■■■ {01,02m} 



{02m, 0l} ■ ■ ■ {02m, 02m} 



Proof. We prove this theorem by induction with m. 

It is true for m = 1 and suppose that {/, f?}^'^-' = {/, 5'}^'^'* for some > 1, we 



shall prove that {f,g} 



(k+i) 



{f,9} 



{2k+2) 
D 



. Let denote a 



;i ■ ■ ■ 



b2k, because 



of (13.71) in the proposition [3] we have: 



{f,9} 



(2fc+2) _ i^[Q02fc+102fc+2/, Q02fc+102fc+25'] 



D 



F[a02fc+l02A:+2, a02fe+102fc+2] 



(3.17) 



Multiplying r.h.s. of fIXTTD by 1 = f|ff}, using the Tanner identities (K2^ . 
(13. 3p and knowing determinant of a skew-symmetric matrix of odd size to be 
zero, F[a(j)2k+i,a(j)2k+i] = 0, we get the r.h.s of (13.171) 



-^[«02fc+l, gg] F[a(j)2k+2f, Q02fc+102fc+2] - i^[«02fc+l , Q02fc+2] -^[Q02fc+2/, a02fc+lg] 

-F[a(p2k+i,0!(l)2k+2] F[a(l)2k+2,a(p2k+i] 

Again, multiplying by 1 = -fjf^, using the Tanner identities (13.21) . the vanish- 
ing determinant of a skew-symmetric matrix of odd size, i.e. F[a(j)2k+2, tt02fc+2] = 
0, and the recursive assumption {M,f}^'^^ = {■u,t>}^''^ we obtain: 



1^ ^|(2A:+2) _ F[af, ag] ^ F[af, a(?!>2fc+i]-F[Q:02fc+2, ag] - F[af, a(/)2fc+2]-F[a02fc+i, aff] 
^ F[a,a] F[a,a]F[acl}2k+i,o:(j)2k+2] 

^ 1^ ^|(2fc) _^ {/,02fc+l}D^^{02fc+2,g}^^^ - {/,02fc+2}^^^{02fc+l,g}^^^ 

{02fc+l, 02fc+2}D^^ 

^ g^ik) ^ {/,02fc+l}^'H02fc+2,ff}^'^ - {/,02fc+2}^'H02fc+l,g}^'^ 

{02fe+l,02fc+2}(^) 

It implies that r.h.s. of eq. (I3.17P is equal {f,g}^'''^^^ which ends the proof. 4 
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Theorems [T] and [2] are main results of this article. 

One may use theorem [2] in proving Jacobi identity and some other algebraic 
properties for Dirac bracket. For example, one can prove the following 

Proposition 5 Suppose (JF, {-, ■}) be skew-symmetric algebra and {(pk, k = 
1, . . . , 2n} be a set of elements from JF such that W{(j)i, . . . , 02n) is invertible. 
Then V/, g & J-": 



^ \W{(f)l,...,(f)2nJ,g)\ ^ F[(j)i ■ ■ ■ (j)2nfg, 01 " " " 02n/5'] 



Proof. Let a = 4>i4>2 ■ ■ ■ 4>2n- From the identity (13. 2p we have 

F[a,a\F[afg,afg] = F[af,af]F[ag,ag] - F[af,ag]F[ag,af] = {F[af,ag])^. 

(3.19) 

Dividing both sides of fl3.19p by {F[a, a]Y (i.e. |iy(0i, . . . , 02n)P ) we obtain 

F[a,a\ XJ^ySD- 



3.3 Jacobi identity 



In |2] , Dirac was struggling to prove the Jacobi identity for his bracket formula. 
He wrote: "I think there ought to be some neat way of proving it, but I haven't 
been able to find it". The Proposition [6] below contains what we believe is just 
that kind of a proof. 

Proposition 6 Let (JF, ■) be a commutative algebra with Lie or Poisson 
bracket {■, ■}. Suppose {0^, k = 1, . . . , 2n} be a set of elements from such 
that ({0j,0j}) is invertible. Then {■,-}d with respect to {(pkYkLi is a Lie or 
Poisson bracket, respectively. 

Proof. Only the Jacobi identity is difficult to verify. Using the theorem [2] and 
the induction principle, it is enough to show that {■, -j*-^^ satisfies the Jacobi 
identity. In order to check the Jacobi identity for {■, ■j'^^^ it is convenient to 
introduce the following symbols: Ai = {/, 0j} , Bi = {g, (pi} , Ci = {h, 0j} with 
i = 1,2 and 0i2 = {0i,02}- Since the Jacobi identity holds for {■, ■} all the 
following sums vanish 



li = {Aug} + {/, Bi} + {^i, {/, g}}, Ji = {A,, h} + {/, d} + {^i, {/, h}}, 
Ki = {Ci,g} + {h, Bi} + {<Pi, {h, g}}, D = {02, ^i} + {A2, 0i} + {/, 0i2}, 
E = {<^2, B,} + {B2, 0i} + {g, 012}, F = {cP2, Ci} + {C2, 0i} + {h, 012}. 
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(3.20) 



Full expansion of Jacohi = {f,{g,h}D}D + {fi', {/i, /}d}d + {h,{f,g}D}D 
produces 39 non-vanishing terms that can be grouped in a polynomial of the 
variable z = {(pu)^^ as follows: 



Jacobi = [{/, {g, h}} + {g, {h, /}} + {h, {/, g}}] + 

[(A2K1 - A1K2) + iB2Ji - B1J2) + (Cih - C2I1)] z + 

[{AiB2 - A2Bi)F + {C1A2 - C2Ai)E + {B1C2 - B2Ci)D] . (3.21) 

Clearly, r.h.s. of (13.211) is equal zero since all its coefficients are zero according 
to ^M). 



4 Applications 



One important class of constrained dynamical systems is characterized by 
K holonomic constraints 0i(g) = 0, where i = 1,---,K. These constraints 
represent a subclass of time-independent constraints 4>i{q,p) = considered 
in this article. In the Dirac approach, these dynamical systems are described 
by a system of 2K constraints 4>i{q) = and (pi{q,p) = {0i,7i} = 0. 

For holonomic constraints, it is convenient to introduce two K x K matrices: 
symmetric S = (Sij) with Sij = and skew-symmetric A = (Aij) with 

Aij = {(pi, (pj}. The matrix W and its inverse C can then be written as 



W 



s 
-5^ A 



and C = 



-1 



S-^AS-^ -S-^ 

5-1 



(4.1) 



In order to compute C one has to invert one symmetric K x K matrix and do 
matrix multiplications twice. Symbolic computation is costly, but numerical 
computation requires only ~ flops (floating-point operations). 

Consider now a constrained model with damping force proportional to the 
generalized velocity. Such a case is described by a metriplectic structure: 



{xi,Xj} = = {Pi,Pj}, {xi,Pj} = 6ij, 
< Xi, Xj >= 0, < pi,pj >= 6ij\i{q,p), where > 0. 

The dissipative constraint matrix = (w^j^, where 
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is a symmetric K x K matrix, and let denote its inverse matrix by 
(W^)~^. The metriplectic Dirac equations for the dynamics governed by / 
{f,'H}D— < /, >D: take the form: 



dpi 

dn 

dq, 

-A,: 



dp, 



■4>k, 



<li = —- Z^iS )jk- 
j,k=l 

pi=-^-f:is-')jk^{4'k,H}+ 

j,k=l 



K 



on sr^ u<4Jj s-^ 



dn 



1=1 



dpi dpi 



4>k 

(4.2) 



Recursive symbohc evaluation of explicit equations for a system having 2K 
constraints is realized by K steps. In each step we deal with only two con- 
straints, e.g (pi and 0j in the i-th step. In order to calculate 2n explicit equa- 
tions of motion subject to 2K constraints, i.e. {a;^, Ti}*^-^) and {pi,'HY^\ we 
have to compute (6n + 3) brackets determined in (X — l)-th step: {xj, 1-LY^~^\ 

{p^,nY''~'\ {xM^''-'\ {x.M^^-^), {piAKY^-'\ {<t>K,nY^-'\ 

{0;., ^F"^) and {4>kAkY''-'^- 

We illustrate our procedure on the model of chain molecule often studied in 
polymer and proteins physics, paying particular attention to the implementa- 
tion of the code for Dirac brackets in symbolic computer algebra system. 

A chain molecules is a constrained system consisting of massive points (or 
spherical balls) attached by rigid massless bonds having fixed length, in d-dim 
space. We are interested in the cases when d — 2 (planar) or 3. The molecules 
interact with each other through a pair potential which depends only on the 
distance between molecules, e.g the Coulomb interaction and/or Lennard- 



Jonnes potential Va — a— -|- e 



(^) -(^) 
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and with an external field 



U{fi). In a real application such a chain is immersed into a fluid matrix, thus 
each of its molecules is subject to an additional frictional force. 

We denote the position of the i-th molecule as and its momentum as pi. 
We will lump all the positions into one vector r = (rl, • ■ ■ , rV) and similarly 
p = (pi, ■ ■ ■ ,Pn)- It is convenient also to use the following notation: the relative 



position of i-th and j-th molecule 



'J' 



the relative position of two 



consecutive molecules (or shortly link vector) Afi 



rj_i_i, the relative 



velocity of two consecutive molecules Avi — £^ — pw, g^j^^j -j-j^g ^^j-j- vector of 



the link vector 



|Ar-| 



. The Hamiltonian for our model reads then 
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mi 



•m2 



mN-i 



Fig. 1. Linear polymer consists of N molecules interacting each other. 



N 



\pi 



2mj 



+ Uin) 



N 



N 



+ E ^^.M = ESr + ^(0- (4.3) 



i=l L^""* J j>i+l i=l 

Putting K = [N — 1), the 2K constraints follow: 



i,(f) = -(|Af,p-/D 



(f, p) = Avk ■ Ark = . 



(4.4) 



Using this notation we can easily evaluate matrix coefficients for all the matri- 
ces in Eq. (14.21) . We found it convenient to collect them in the tabled! where 



the 6j, Cj, Oj, b\^^ and cl^' (for isotropic friction Xi(^d-i)+i = ■ ■ 
which is the frictional coefficient for i-th molecule) are given as 



A,- 



i{d-l)+d 



A,; 



cos(aj) , where cos(ai) = — • e^+i 



(mj + mi+i) 2 
Ci = An 



1 1 

+ 



rrii mj+i 



t2 

H 1 O-i 



Afi ■ Avi+i - Avi ■ Afj+i 



m 



^^iiAfj • Afi+i = — 



Ai+i 



m 



kk+i cos(Qj) 



Thus, the matrices S*, S^^^ are symmetric tridiagonal, while A is skew- 
symmetric tridiagonal, shown in the table [2l For homogeneous polymer in ho- 
mogeneous environment, consisting of identical molecules, U = I and rrii = m, 
all formulas on elements of S, S^^^ become even simpler: 



— , bi = — cos(aj), and c,- = — ^, b) = — ^ 
m m 



cos a,- . 



(4.5) 
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Condition 


Sij = {4>i,4>j} 


Aj = Wi,4>j} 


sip =<k4>j> 


> 1 











j = i+l 


bi 


«i 




j = i 


Ci 







j = i-l 









Table 1 

Elements of the matrices S, A and S"^^^ 



S 



ci 6i ••• 

6i C2 62 '■• : 

■•• ■•• ■•• 

: ■•• bx-i 

••• bx^i CK 



and A 



ai ••• 
-ai 02 '■• : 

■•• ■•• ■•• 

: ■•• ax-i 

•••0 -ax-i 



Table 2 

Symmetric and skew-symmetric Tridiagonal Matrices S and A 

Though the tridiagonal matrices have been considered numerically for years, 
the explicit analytic formulas for elements of the inverse matrix of a tridiagonal 
matrix are known only in some special cases [18J: hi = h and Cj = c. Here we 
propose a general expression for elements of S^^. Details of the derivation of 
that formula are given in the Appendix A. 

Let S{1, — 1) be the top left {i — 1) x [i — 1) matrix containing rows 
and columns {1, . . . , i — 1} of S and S{j + 1, ■ ■ ■ , K) be the bottom right 
{K — j) X [K — j) matrix containing rows and columns {j + 1, . . . , K} of S, 
we get the following recursive formula: 



for i < j, and is symmetric. Since both matrices S and S"^^^ have a similar 
form, we can use the formula (14. 6p in calculating their inverse. 

Furthermore, for K > n > I > 1, the \S{1, ■ ■ ■ ,n)\ is calculated from the 
recursive relation: |S'(0)| = 1, |5'(/)| = q, \S{1, ■ ■ ■ ,n)\ = c„|S'(/, ■ ■ ■ , n — 1)| — 
6t,|5(/,---,n-2)|. 
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With the formula (14.60 , it is easy to show that the inverse matrix of a symmet- 
ric tridiagonal matrix is one-pair matrix. Numerically it can be computed fast 
with 0{N) complexity cost, and with modest memory usage. Since the recur- 
sion relation fl4.6p is rather involved, we can only calculate the Dirac equations 
via recursion. More technical details are presented in our paper posted on the 
arxiv page. 

Discussion 

We have implemented our formalism using the package Mathematica version 
5.2 and 6.0, the computer algebra system, both for symbolic and numerical 
calculation, and measured the CPU time needed in computing explicit analyt- 
ical r.h.s. of (14.21) in two ways: one based on the formula (14. 6 p and the other 
based on the recursion relation (I3.15p . All computation have been done on an 
ordinary PC (with dual core processor 1.6 GHz and 1GB RAM) running MS 
Windows XP and Linux FC6. The symbolic computing time for one pair of 



Time (in Sec) 




Constraints 



4 3 12 16 

equations in 3-dim, after using least square interpolation, seems to grow with 
the number of constraints like 0.028 e°'^^^ and as 0.00046 e^ ''^'^ for method 
inverting triangular matrices and using recursive formula, respectively. Con- 
sequently, the recursive formula is reasonably good only for systems with less 
than 12 constraints. Since the computing time in both methods grow expo- 
nentially in the number of constraints, computing explicit analytical Dirac 
equations seems to be inapplicable for very long chains. However, fast algo- 
rithm for numerical inversion of tridiagonal matrices does exist and has a 
complexity 0{N). Thus, Dirac finite difference equations for long chains are 
computable. 

Having explicit equations of motion one can solve them numerically either by 
standard explicit /implicit Runger-Kutta algorithm or standard Mathematica's 
ODE solver NDSolve. 
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Another important issue is that alternatively to the system of equations (14. 2p . 
one can consider the following system: 



Pi 



dn 

_m 

dqi 



K 

j,fc=i 



•jk^{(i>k,'H} 
dqi 



dpi 



K 



j,k=l 



(4.7) 

dpi dpi 



Since constraints are Casimir elements regarding to Dirac bracket, any solution 
of (14. 2 p with initial conditions satisfying all constraints, automatically satisfies 
all constraints for all time. Therefore it must also be a solution of (14.71) . 

This fact and the uniqueness of solution (locally) imply that two systems 
(14. 2 p and (14.71) are equivalent. In our tests, symbolic computation for the 
latter is 6-7 times faster than for the former. Moreover, for non-dissipative 
mechanical systems, the latter is exactly the system of equations obtained from 
the Lagrange Multiplier Method (LMM), eq. flA.Sp in the Appendix |X1 Though 
these two systems are mathematically equivalent, they are not equivalent for 
numerical algorithms approximating solution, which means that errors grow 
differently for each of them even if using a common numerical algorithm. Errors 
in computing approximate solution of the LMM-like eq. (14.70 or (lA.Sp . always 
grow faster than those of the Dirac-like eq. (14. 2p . We studied numerically the 
violation of energy and bond length constraints for a particular polymer with 
one fixed end, eg. A^-pendulum described in the Appendix |Bl These numerical 
results are presented briefly in the figure [2l In summation, standard numerical 




(a) Energy cal- (b) Energy (c) Sum of con- (d) Sum of 
culated from eq. calculated straints errors constraints 
(|4.2p and ()4.7p from (jA.SP calculated from errors calcu- 

eq. (j4.2p lated from eq. 

Fig. 2. Numerical test: Energy and Constraints errors for 4-pendulums dynamics 
described by the Hamilton-Dirac eq. (j4.2p . simplified Dirac (14. 7p and LMM (jA.5P 
using default numerical algorithm NDSolve. For simplicity we have chosen a system 
consisting of 4 equal masses which are in the axis x at the beginning, and whose 
initial velocities have random values satisfying constraints' equations. 

algorithms seem to work well with Dirac-like equations. To deal numerically 
with LMM-like equations, we recommend to use either constrained algorithms 
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(eg. SHAKE, LINGS) or other advanced symplectic/poisson ones, which have 
been developed recently. 

Although in the simulation, polymers with nearly constant bond length, called 
stiff bead-spring chains, are more often considered than those with rigid con- 
stant length, named bead-rod chains, the matrix S which has been carefully 
studied here, is closely related to the metric potential U = ^kTlog{\S\) in the 
statistical mechanics of Polymers |23j . 

The application of bracket formalism to the non-linear many particle models 
is possible by time consuming. We have looked at the possibility of using our 
method to obtain a set of analytical equations and simulate mechanics of the 
caricatured human body [19]. 

Instead of models for body dynamics such as inverted pendulum [20] , or elastic 
string [2T] are used, we used skeletal humanoid consisting of 13 material points, 
fig. [3l We found that symbolic calculation each pair of explicit analytical 
equations for humanoid takes app. 9 minutes using formula (14.61) for inverting 
matrix S, of uninterrupted Mathematica performance in PC. 



Fig. 3. Humanoid is a (dissipative) constrained dynamical system with 24 phase 
space constraints. This is an example of non-Hnear chain. 



5 Conclusions 

In this article we have reviewed a geometric construction of Dirac-like brackets 
and proved recursive character of such brackets. We showed that computing 
explicit dynamical equations based on these brackets may be difficult, but 
it is possible to produce analytical equations even for systems with many 
constraints. 

We have applied here the Dirac procedure for metriplectic mechanical models 
with finite degrees of freedom, but in our previous work we have shown its use- 
fulness for continuous models [14j , for example incompressible hydrodynamics 
[22] • Fixman [23] have used constraints approach in formulation of statistical 
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mechanics of various polymer models. The fact that constraints can then be 
visualized as a kind of temperature dependent potential is not unusual. Fix- 
man and others have restricted their procedure to the equilibrium calculations. 
Our formalism allows us to go beyond the equilibrium application and see the 
form of the constrained Liouville equations, modifications in the dynamical 
modes coupling due to the presence of constraints and possible the role of the 
constraints play in removing the singularities appearing in low dimensional 
systems statistical mechanics. For example, the fact that the transport co- 
efficients, like viscosity, thermal conductivity and diffusion coefficient do not 
exists in d = 2, can be modified by presence of the constraints in a fashion 
analogous to that mentioned in |24j . 
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A Lagrange Multiplier Method 

The purpose of this section is to show that computing explicit analytical equa- 
tions in the Lagrangian formalism is equally difficult as in the Dirac formalism. 

For simplicity, suppose that all constraints of the form: (j)k{(l) = 0, k = 
1, . . . ,K and q = (gi, . . . , Lagrangian of constrained system is a sum of 
unconstrained Lagrangian and a linear combination of constraints: C{q, q) = 
^o(?i <i) ~ Y.k=i ^k<Pk{<l)- The Euler-Lagrange equations read 

d_ fdC\ _dC_^ 
dt \dq J dq 

Suppose Lagrangian of the form Cq = T{q) — V{q) = ^q'^ M q — V{q), with in- 
troducing conservative force F = — the Euler-Lagrange equations become 



Mq = F-J2Xk^ = F-BX, (Al) 

fc=l 



21 



where B = (Bik) is a nxK matrix whose elements Bik = Since = 0, 
all first and second time derivatives of (pk vanish: 



where Gk = YA,j=i dq'dq^ '^^r Substituting for q = M'-^lF — BX], derived from 
(EH), in (D) we get: ' 



= G + B^M-^[F - BX], (A.4) 

here G = (Gk), X = (Afc) are column vectors K x 1 and F = (Fj) is a 
column vector n x 1. Therefore, [G + B'^Ar^F] = {B^M-^B)X or A = 
{B^M-^B)-^[G + B^M-^F]. Substituting this back to we get exphcit 

constrained equations: 



Mq = F - B{B^M-^B)-^[G + B^M-^F] . (A.5) 

Thus, for achieving explicit equations in the Lagrangian formalism, it is also 
necessary to compute analytical inversion of the K X K matrix (B^M'^B) 
which is exactly equal the matrix S in the Dirac approach where the Hamil- 
tonian obtained from the Legendre transformation: 7i = pq — C with P = 



B N-pendulum in d dimensional space 



We denote the position of the i-th mass as rj = {xci{i~i)+i, ■ ■ ■ ,Xdi), its mo- 
mentum as Pi = {pci{i-i)+i, . . . ,Pdi), the relative position of i-th and j-th 
mass fij = fi — fj, the relative position of two consecutive masses (or shortly 
link vector) Afi = fi — "Ti+i, the relative velocity of two consecutive masses 
Avi = ^ — Pi+i- and the unit vector of the link vector Cj = i^^. 



B. 1 Hamilton- Dirac description for N-pendulum 



The Hamiltonian is given by T-L{f,p) = J2i=i + Xd. 
class constraints follow: 



^ + gmi Xdi , and 2N second- 
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'N-1 



Fig. B.l. N-pendulum is a constrained system with length constraints, which can 
be viewed as a hnear polymer with fixed end. 



i(E?=ix|-/i2) = l(|riP-/2) forfc = l 



4>k{r,P) = {0fc,^} = 



for 1 < A; < 



ri ■ Vi for /c = 1 , 
Avk ■ Affc for 1 < A; < AT . 



(B.l) 
(B.2) 



B.2 Lagrange Multiplier Method for N-pendulum 



The Lagrangian is given by C{f,p) = J2iLi 
constraints follow: 



2m,- 



g rrii Xdi 



and A^ length- 



^ ; \[T.U^]-ll] = \{\ri?-ll) forfc = l 



i(|Affc|2-/,2) 



for 1 < A; < A^ . 



(B.3) 



In order to calculate explicit eq. ( 1A.5I) we need to calculate explicit elements 
of where S follows: 



b'^m-^b = S 



ci 6i 

bi C2 62 : 

62 C3 63 ■•• : 

: ■•. ■•. ■•. ■•. 

: • • ■ • CN-i bN-i 

bN-i cn 



(B.4) 



here 
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1 -Ar 



hh 



mi mi 
Ari-i-Ari _ kk+i 



cos(ai) for z = 1 , 

cos(aj) for 1 < i < — 1 



mi 



for ^ = 1 , 
f ^ + If for 1< z < AT . 



(B.5) 



C Symbolic Inversion of Symmetric Tridiagonal Matrices 



In this section we discuss problem of symbolic inversion general symmetric 
tridiagonal matrix whose explicit form is given in flC.l|) . 



S 



(C.l) 



ci 6i 

bi C2 62 : 

62 C3 63 ■•• : 

: ■•. ■•. ••• ■•. 

: ■ • ■ • ck-1 bx-i 

CK 



Notation: Let M — {Mij) be a matrix. Define as M[ii, . . . ,ip;ji, . . . ,jq ) the 
matrix consisting of elements Mj j where i G {ii, . . . ,ip} and j G {ji, . . . ,jg}. 
In the case {ii, . . . ,ip} = {ji, . . . ,jg} instead writing M{ii, . . . ,ip;ii, . . . , ip) 
we will write M(zi, . . . ,ip). 

Definition 3 A n x n symmetric matrix Q is called an one-pair matrix if 
its elements are products of components of two vectors u = {ui, . . . ,Un) and 
w = (wi, . . .,Wn), i.e. 



Q 



UiWj for i < j 
UjWi for i > j . 



(C.2) 



C.l Direct computation 



Since {S = (— 1)*+-'|S'*^'''*^|/ det S, in order to compute elements of S* ^ one 
has to compute the determinant 15*1 = detS* and the co-factor {—iy~^^S^^'^^\. 
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First, if denote the determinant of k x k symmetric tri diagonal matrix by 5*^, 
then Sk can be calculated recursively as follows: 



5*0 — 1, 5*1 — Ci, Sk — CkSk^i — b1_^Sk-2- (C-3) 

Second, since S is tridiagonal, the S^^'^^ has three decoupled sub-blocks on the 
main diagonal with the order {i — 1), \j — i\ and {K — j). 

Denote the determinant of a matrix (Sab) whose indexes a, b belong to the set 
{ii,Z2, ■ ■ - ^ip} by \S{ii,i2, ■ ■ ■ ,ip)\, with i < j we have 



For K > k > I > 1, the \S{1, . . . , k)\ is computed from recursive relation: 
\Sm\ = 1, \S{1)\ = Q, \S{1, . . . , fc)| = Ck\S{l, . . . ,k~l)\-bU\S{l, ...,k-2)\. 

Introducing df^ = \S{1, ...,/ + i — 1) | with I + i — 1 < K, the sequence df^ for 
fixed value of / is determined by the recursion: 



4'^ = 0, d? = ci, 4?i = - • (C.5) 

The elements of the inverse matrix are calculated as follows: 

{S'\, = {-iy+^ bA+i ■ ■ -fe,-! , for ^ < J . (C.6) 

If bi 7^ 0, one can introduce 



= -^15(1, . . . , ^ - 1) I bA+i ■ ■ ■ bK-i (C.7) 

and express elements of S^^ by two vectors (-Ufc) and (wk) 

{S^^)ij = UiWj for i < j . (C.9) 

Thus, we have proved that the inverse of a symmetric tridiagonal matrix is 
one-pair matrix. The reverse statement remains true. Eq. (IC.SP and (1C.6P seem 



to define the most effective algorithm for computing elements of 5* ^. 
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C.2 Block diagonalization 



This method based on the observation that for a symmetric tridiagonal matrix 
S it is easy to find a sequence of upper triangular (non- symmetric) tridiagonal 
matrices Uk, k = 1, - ■ ■ , K—1, with the main diagonal {1, . . . , 1, Xk, Zk,l, ■ ■ ■ , 1} 
and its upper neighbour diagonal {0, . . . , 0, y/j, 0, . . . , 0}, i.e. 



Uk- 



1 

■•. 







Xk Vk 
Zk 



'■■ 
1 



(C.IO) 



such that VA;, 1 < k < K - 1: 



{Ui... Ukf S{Ui... Uk) 



4+1 : 

Pk+l ' ■ '■ 

(3k+i Cfe+2 bk+2 '■ ■ 

: ■■• bk+2 '■■ '■■ 

: ■•• ■•. ■•. bx-i 

•••0 bx-i CK 



(C.ll) 



Thus, their product U — U\ - ■ ■ Uk-i is the upper triangular satisfying: 
S U = I. This implies that the inverse matrix of S" is a product of U 
and its transposition, — UU^ . With convention zq — 1, elements of U 
follow: 



Ui,j = < 



Xj{yi ■ ■ ■ yj_i)zi_i for i < j - 1 

XjZj^i for i = j 
for i > j . 



(C.12) 



In order to calculate Uk one has to solve recursively a system of 3 quadratic 
equations 
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Xk 




O-k Pk 




Xk Vk 


Vk Zk 




Pk dk 




Zk 



(C.13) 



where 



ai = Ci, Pi = bi, di = C2, and for A; > 1: 

O-k = 1, Pk = bkZk-l = bk. 



dk-l 



ak-idk-i - iPk-iY" 



dk — Cfc+i- 



(C.14) 



In each stage, the quadratic system (IC.lSp has 4 solutions but for our purpose 
it is enough to consider only one among them 



1 



Xk 



Vk 



Pk 



Zk 



\J ak[akdk - [PkY] \lakdk-{PkY 
Express the elements of the inverse matrix = UU'^ respecting flC.121) 



(C.15) 



K-1 K-1 

{S )i,j= E ^i,kUj± = E Ui^kUj^k 
k=l fc=max{ij'} 
K-1 

= E {,Xk)'^{yi. . .yk-i){,yj ■ ■ ■yk-l)z^-lZj^l. 

fc=max{i j'} 

Note that in flC.161) . recursion appears only in the expression of Pk'- 



(C.16) 



Pi = hi, P2 = b2. 



Cl 



C1C2 - bi 



,2 ' 



Pk = h 



Ck - {Pk-lf 



for A; > 2. (C.17) 
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